1 Poisson Events

1.1 Libraries

library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
layout(matrix(1:1, nrow=1))

1.2 Parameters and risk

censoredProb <- 0.0002
timeSpan <- 20
timeInterval = 0.01
InitialPopulatoin <- 1000
ContBetaRate_1 <- 0.0002
ContBetaRate_2 <- 0.000001
BinBetaRate_1 <- 0.001
BinBetaRate_2 <- 0.003
BaselineHazard <- ContBetaRate_1
betaRates <- c(BaselineHazard,ContBetaRate_1,ContBetaRate_2,BinBetaRate_1,BinBetaRate_2)
BaseVar <- rep(1,InitialPopulatoin)
ContVar_1 <- runif(InitialPopulatoin)
summary(ContVar_1)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000764 0.2350957 0.4777109 0.4771627 0.7061444 0.9996590
ContVar_2 <- rnorm(InitialPopulatoin,50,10)
ContVar_2[ContVar_2 < 1] <- 1 
summary(ContVar_2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   20.34   44.05   50.62   50.50   57.16   82.34
BinVar_1 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_1)
## BinVar_1
##   0   1 
## 628 372
BinVar_2 <- rbinom(InitialPopulatoin,1,0.35)
table(BinVar_2)
## BinVar_2
##   0   1 
## 661 339
dataFeatures <- as.matrix(cbind(BaseVar,ContVar_1,ContVar_2,BinVar_1,BinVar_2))
hazardRate <- as.numeric(dataFeatures %*% betaRates)
summary(hazardRate)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0002367 0.0003613 0.0013161 0.0017349 0.0033288 0.0044498

1.2.1 Getting the events and time to event

aliveSet <- c(1:InitialPopulatoin)
eventSet <- numeric(InitialPopulatoin)
timetoEvent <- numeric(InitialPopulatoin)

for (time in c(1:(timeSpan/timeInterval)))
{
  randProb <- runif(length(aliveSet))
  Iscensored <- randProb <= censoredProb
  Isevent <- randProb <= (1.0-exp(-hazardRate[aliveSet]))
  Iscensored <- Iscensored & !Isevent
  eventSet[aliveSet] <- Isevent
  timetoEvent[aliveSet] <- time*timeInterval-timeInterval/2
  isCensoredOrEvent <- Iscensored | Isevent
  aliveSet <- aliveSet[!isCensoredOrEvent]
#  cat(length(aliveSet),"(",sum(isCensoredOrEvent),",",sum(Isevent),",",sum(Iscensored),")\n")
}
timetoEvent[aliveSet] <- time*timeInterval + timeInterval/2
summary(timetoEvent)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.005   1.915   5.685   8.504  16.413  20.005
table(eventSet)
## eventSet
##   0   1 
## 221 779
pevent <- (1.0-exp(-hazardRate))
summary(pevent)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0002366 0.0003612 0.0013152 0.0017323 0.0033233 0.0044399
simulatedDataFrame <- as.data.frame(cbind(status=eventSet,time=timetoEvent,pevent=pevent,dataFeatures))

1.3 RRplots()

plotTimeInterval <- 2.0

hazard <- -log(1.0-simulatedDataFrame$pevent)
hboost <- plotTimeInterval/timeInterval
pvalue <- 1.0-exp(-hboost*hazard)


rdata <- cbind(simulatedDataFrame$status,pvalue)
summary(rdata[,2])

Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04623 0.06971 0.23143 0.26261 0.48612 0.58933

table(simulatedDataFrame$status)

0 1 221 779

RRAnalysisCI <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=simulatedDataFrame$time,
                     title="Simulation",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=plotTimeInterval)

par(op)

1.3.1 By Risk Categories

obsexp <- RRAnalysisCI$OERatio$atThrEstimates
pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 779 725.2 835.7 770.5 1.011 0.941 1.08 0.759
low 178 152.8 206.2 167.8 1.061 0.910 1.23 0.418
90% 29 19.4 41.6 28.1 1.032 0.691 1.48 0.850
80% 572 526.1 620.9 575.6 0.994 0.914 1.08 0.900 
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))

plot(obsexp$Expected,obsexp$Observed,
     xlim=c(minx,maxx),
     ylim=c(minx,maxx),
     main="Cal. Expected vs Observed",
     ylab="Observed",
     xlab="Expected",
     col=rainbow(nrow(obsexp)),
     log="xy")

errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

1.3.2 Time to Event Analysis

isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Observed Time",
                   at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames") 
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
  O:Low Risk < 0.080 O:0.080 <= Risk < 0.086
1Q 3.57 3.25
Median 9.38 5.44
3Q 13.34 11.45
Table continues below
  O:High Risk >= 0.086 E:Low Risk < 0.080
1Q 1.13 13.4
Median 2.56 15.1
3Q 5.71 16.6
  E:0.080 <= Risk < 0.086 E:High Risk >= 0.086
1Q 11.4 1.46
Median 11.5 1.52
3Q 11.7 3.66 
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
     xlab="Mean Expected Time",
     ylab="Mean Observed",
     xlim=c(minv,maxv),
     ylim=c(minv,maxv),
     main="Estimated Time to Event",col=rainbow(length(classnames)))

lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)

1.3.3 Risk Calibration

op <- par(no.readonly = TRUE)


crdata <- cbind(simulatedDataFrame$status,pvalue,simulatedDataFrame$time)

#calprob <- CalibrationProbPoissonRisk(crdata,timeInterval=plotTimeInterval)
calprob <- CalibrationProbPoissonRisk(crdata)

( 20.28598 , 11004.11 , 1.935431 , 779 , 1494.978 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
6.6 10.2 20.3

1.3.4 After Calibration

h0 <- calprob$h0

caldata <- cbind(simulatedDataFrame$status,calprob$prob)


rrAnalysisTrain <- RRPlot(caldata,atRate=c(0.90,0.80),
                     timetoEvent=simulatedDataFrame$time,
                     title="Cal. Simulation",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=calprob$timeInterval)

1.3.5 By Risk Categories

obsexp <- rrAnalysisTrain$OERatio$atThrEstimates
pander::pander(obsexp)
  Observed L.CI H.CI Expected O/E Low Upper pvalue
Total 779 725.2 835.7 772.4 1.009 0.939 1.08 0.801
low 178 152.8 206.2 168.2 1.058 0.908 1.23 0.441
90% 29 19.4 41.6 28.2 1.030 0.690 1.48 0.850
80% 572 526.1 620.9 577.0 0.991 0.912 1.08 0.851 
maxx <- 1.1*max(c(obsexp$Observed,obsexp$Expected))
minx <- 0.9*min(c(obsexp$Observed,obsexp$Expected))

plot(obsexp$Expected,obsexp$Observed,
     xlim=c(minx,maxx),
     ylim=c(minx,maxx),
     main="Cal. Expected vs Observed",
     ylab="Observed",
     xlab="Expected",
     col=rainbow(nrow(obsexp)),
     log="xy")

errbar(obsexp$Expected,obsexp$Observed,obsexp$L.CI,obsexp$H.CI,add=TRUE,pch=0,errbar.col=rainbow(nrow(obsexp)),cex=0.75)
lines(x=c(1,maxx),y=c(1,maxx),lty=2)
text(obsexp$Expected,obsexp$Observed,rownames(obsexp),pos=2,cex=0.75)

1.3.6 Time to Event Analysis

isevent <- RRAnalysisCI$timetoEventData$eStatus == 1
exptime <- boxplot(RRAnalysisCI$timetoEventData$expectedTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Expected Time",plot=FALSE)
obstime <- boxplot(RRAnalysisCI$timetoEventData$eTime[isevent]~RRAnalysisCI$timetoEventData$class[isevent],
                   xlab="Class",
                   ylab="Time",
                   main="Observed Time",
                   at=exptime$stats[3,],plot=FALSE)
classnames <- attr(RRAnalysisCI$timetoEventData,"ClassNames") 
timesdata <- cbind(obstime$stats[c(2,3,4),],exptime$stats[c(2,3,4),])
rownames(timesdata) <- c("1Q","Median","3Q")
colnames(timesdata) <- c(paste("O",classnames,sep=":"),paste("E",classnames,sep=":"))
pander::pander(timesdata)
Table continues below
  O:Low Risk < 0.080 O:0.080 <= Risk < 0.086
1Q 3.57 3.25
Median 9.38 5.44
3Q 13.34 11.45
Table continues below
  O:High Risk >= 0.086 E:Low Risk < 0.080
1Q 1.13 13.4
Median 2.56 15.1
3Q 5.71 16.6
  E:0.080 <= Risk < 0.086 E:High Risk >= 0.086
1Q 11.4 1.46
Median 11.5 1.52
3Q 11.7 3.66 
minv <- min(c(exptime$stats[2,],obstime$stats[2,]))
if (minv<0.001) minv <-0.001
maxv <- max(c(exptime$stats[4,],obstime$stats[4,]))
errbar(exptime$stats[3,],obstime$stats[3,],obstime$stats[2,],obstime$stats[4,],log="xy",
     xlab="Mean Expected Time",
     ylab="Mean Observed",
     xlim=c(minv,maxv),
     ylim=c(minv,maxv),
     main="Estimated Time to Event",col=rainbow(length(classnames)))

lines(x=c(minv,maxv),y=c(minv,maxv),col="black",lty=2)
legend("topleft",legend=classnames,lty=c(0),pch=c(16),col=rainbow(length(classnames)),cex=0.80)